3D simulations of linearized scalar fields in Kerr spacetime 
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We investigate the behavior of a dynamical scalar field on a fixed Kerr background in Kerr- 
Schild coordinates using a 3+1 dimensional spectral evolution code, and we measure the power- 
law tail decay that occurs at late times. We compare evolutions of initial data proportional to 
/ (r)Y( m (0 , <f>) , where Yi m is a spherical harmonic and (r,9,<f>) are Kerr-Schild coordinates, to that 
of initial data proportional to /(r B L)^m(#BL, 4>)> where (Ybl, #bl) are Boyer-Lindquist coordinates. 
We find that although these two cases are initially almost identical, the evolution can be quite 
different at intermediate times; however, at late times the power-law decay rates are equal. 
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I. INTRODUCTION 

The propagation of classical scalar fields in a fixed 
black hole spacetime has been studied extensively ever 
since the work of Price Q that described the behavior 
of such fields in the Schwarzschild geometry. Higher-spin 
fields, such as linearized gravitational perturbations, be- 
have qualitatively similar to the zero-spin case, and there- 
fore scalar fields are often used to gain insight into more 
general situations. Although the behavior of scalar fields 
in Schwarzschild spacetime is well-understood, the situa- 
tion for a Kerr background geometry is still under active 
investigation and has been a topic of some controversy 
(see, e.g., 001 and references therein). 

The evolution of a scalar field in curved spacetime is 
governed by the massless Klein-Gordon equation 

Of>(x,y,z,t)=0, (1.1) 

where ib is the value of the scalar field and □ is the 
dAlembertian operator in curved spacetime. Accord- 
ing to no-hair theorems, the only nonsingular time- 
independent solution to Eq. (|l.l(l in a black hole back- 
ground is tp — everywhere, and furthermore, if tp ini- 
tially varies in time or space it will evolve until it reaches 
this time- independent solution 0. When observed at a 
fixed spatial location as a function of time, the evolu- 
tion of a scalar field in a black hole spacetime consists of 
three distinct phases, as shown in Fig.^ The first stage 
is the initial burst, which is determined by the initial 
conditions imposed on the scalar field. The second stage 
is the quasinormal ringing phase, during which outgo- 
ing waves interfere with incoming waves that backscatter 
off the black hole's potential well. During this phase ib 
oscillates and decays exponentially, and can be written 
as a sum of terms of the form e luJnt for a discrete set 
of complex eigenfrequencies uj n . During the third stage, 
or tail phase, ib depends on incoming radiation that has 
been backscattered off the spacetime curvature at large 
distances. During the tail phase, the scalar field decays 




t/M 

FIG. 1: The evolution of a scalar field with an initial Yio 
angular dependence in Schwarzschild spacetime. Plotted are 
the L2 norms of i/j and t/iona spherical surface of fixed radius 
r, defined by (||/||i2) 2 = (l/4.7r) J f 2 dQ.. Here the integration 
is taken over a surface at r = 11. 9M. The duration of the 
initial burst is about 50Af. After the initial burst, the scalar 
field settles into the quasinormal ringing phase until about 
200Af when the tail phase begins. 



as a power law, yb oc , for the case of a Schwarzschild 
background, and there is good analytical and numer- 
ical [3, Q evidence that it decays as a power law for a 
Kerr background as well. 

For a scalar field in Schwarzschild spacetime, the 
power-law decay in the tail phase is computed using the 
spherical harmonic decomposition of the initial data. The 
amplitude of each Yi m mode present in the initial data 
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will eventually decay like t~^ 2e+3) at late times dHHH, 
assuming that ip initially falls off quickly enough at in- 
finity 1 . If one measures ip at a single point in space or by 
some other method that does not select specific spher- 
ical harmonic components, the decay rate measured at 
late times will be determined by the smallest £ present in 
the initial data, because this is the most slowly decaying 
mode. 

The late-time behavior of the scalar field becomes more 
complicated in Kerr spacetime because of the lack of 
spherical symmetry. Although axisymmetry prevents 
mixing of spherical harmonics with different to values, 
harmonics with different values of £ no longer evolve in- 
dependently. Because of this mode mixing, if the ini- 
tial data are proportional to a pure spherical harmonic 
Ye m i the- evolution should produce spherical harmon- 
ics with different values of £, and in particular smaller 
values of I. It is not unreasonable to assume the same 
power-law time dependence as the Schwarzschild case, 
namely ip ~ t~^ 2l+3 \ because the tails are due to radia- 
tion backscattered off the weak-field asymptotic region of 
spacetime. Given this assumption, the late-time behav- 
ior of the scalar field should be dominated by the smallest 
value of £ that is produced by mode mixing, since this is 
the most slowly decaying mode. Given that £ > m and 
that parity is conserved (i.e. the equatorial symmetry 
of the initial data is preserved), the lowest-order spher- 
ical harmonic that may be generated from initial data 
proportional to Yi omo is £ = too if £q — mo is even and 
£ = too + 1 if £q — to is odd 2 . Therefore, according to 
this simple picture one would expect at late times 



( y f -(2m +3) 
^ I Y [ma + 1)mo t -Wm 0+ X )+ 3) > lQ _ mQ qM CM 

However, analytical work by Hod [Tol | predicts different 
behavior for scalar fields in Kerr spacetime. According to 
Hod's analysis, the late-time decay rate does not just de- 
pend on the lowest multipole index £ permitted by parity 
and axisymmetry; the initial value £q also plays a role: 



ip oc 



Ye +-( 2f o+3) 

V +-(^o+mo + l) 

1 rngmo L 



£ - to < 2, 
£q — too > 2 (even) , 
Y {ma+ i )ma r< Wm »+ 2 ), £ - m > 2 (odd). 



(1-3) 

This is a deeply surprising result, for it implies that the 
generated modes somehow "remember" the properties of 
the initial data that created them. 

It has been recently argued by Poisson Q that both 
Eq. 1|1.3[) and the simple picture leading to Eq. 1|1.2[) are 



1 If ip at large distances is initially a static solution of Eq. 11.11 . 
then the late-time decay rate lj is t~ ( 2 *+ 2 ). Note that all static 
solutions of Eq. 11.11 that are regular at infinity diverge at the 
horizon 00 

2 The amplitude of this lowest-order mode may, however, turn out 
to be zero, in which case the decay rate would be determined by 
the lowest-order mode with nonzero amplitude. 



valid descriptions of the late-time dependence of scalar 
fields in a rotating spacetime that is weakly curved ev- 
erywhere; the difference is merely the choice of spatial 
coordinates (r, 9, </>) used when setting the initial data. 
Poisson assumes a metric equal to Minkowski space in 
spherical coordinates (r, 9, <p) plus a stationary nonspher- 
ical perturbation which he treats to linear order. For ini- 
tial data proportional to f(r)Yi omo (9, <p) for some func- 
tion /(r), he finds that the field decays like t~( 2£ o+3) 
(there is no mode mixing to first order in the pertur- 
bation). He then repeats the calculation using initial 
data proportional to f(r')Y£ omo (9',(j)'), where (r',6',(f)') 
are spheroidal coordinates defined by 

r 2 sin 2 9 r 2 cos 2 9 
■> — + — = !. (1-4) 



r' cos ( 



= rcos6>, (1-5) 

<t>' = 4, (1.6) 

for some constant a. In this case he finds that the modes 
mix because of the nonspherical coordinates, and the 
scalar field decays according to Eq. I|1.3|) . 

Poisson then argues that since radiative falloff is es- 
sentially a weak-field phenomenon, similar conclusions 
should be true for scalar fields in a Kerr background, 
so that coordinate effects would account for the discrep- 
ancy between Eq. Ijl.2|l and Eq. ljl.3|l . Consider initial 
data proportional to j '{r)Y^ oma {9 , <p) where (r,9,(p) are 
any coordinates in which the weak-field limit of the Kerr 
metric is spatially isotropic. Then the only mode mixing 
will be due to the strong-gravity region at early times, 
and at late times, when the scalar field probes only the 
weak-gravity region, each mode that was generated by 
the mixing will decay like £~( 2£ + 3 ). One therefore ex- 
pects Eq. I|1.2|) to hold. Now consider initial data propor- 
tional to f(r')Ye omo (9' , <j)') where (r', 9', <f)') are any coor- 
dinates in which the weak-field limit of the Kerr metric is 
spheroidal. Such coordinates include Boyer-Lindquist co- 
ordinates, the coordinates used in Hod's analysis, which 
in the weak-field limit reduce to flat space in spheroidal 
coordinates (/, 9' , </>') with the parameter a in Eq. I|1.4|) 
equal to the Kerr spin parameter. For such initial data, 
if strong-gravity mode mixing can be ignored relative to 
the mode mixing resulting from the spheroidal coordinate 
system (this key assumption is discussed in more detail 
in Section TlV CI) . then the scalar field should behave ac- 
cording to Eq. HI. 3D . The "memory" effect implied by 
Eq. H1.3[) . according to this argument, is due to coordi- 
nates and not physics. 

Surprising and seemingly contradictory results have re- 
sulted not only from analytic studies of this problem but 
also from numerical simulations. Early simulations 
considered cases for which Eq. I|1.2|l and Eq. (|1.3|l agree, 
and were consistent with both predictions, but a more 
recent simulation yielded the puzzling result that a 
scalar field initially proportional to Y40, the lowest mul- 
tipole mode for which Eq. (|1.2|l and Eq. I|1.3|l differ, de- 
cays approximately like i -5 5 , in conflict with both pre- 
dictions. Most recently, a 2+1 simulation of an initial 
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Y40 mode using ingoing Kerr coordinates agrees with 
Eq. i|1.2f> to high accuracy. 

Here we solve the scalar wave equation in a fixed Kerr 
background in Kerr-Schild coordinates using a 3+1 nu- 
merical evolution code. We reproduce the known funda- 
mental quasinormal frequency and the known tail falloff 
behavior for the Schwarzschild case. We find that for a 
black hole with nonzero spin, when we choose initial data 
proportional to f{r)Y^ (9, </>), where (r, 0) are Kerr- 
Schild coordinates, we find late-time tail behavior con- 
sistent with Eq. HI. 21) . 

We then choose a different set of initial data pro- 
portional to /(>bl)1«)(#bl, 4>), where (Y B l,#bl) are 
the Boyer-Lindquist coordinates. Note that the Boyer- 
Lindquist coordinates are spheroidal in the sense dis- 
cussed earlier and the Kerr-Schild coordinates are spher- 
ical, and that the transformation between (r, 9) and 
(?"bl,#bl) is the same transformation i|1.4H1.5fl used by 
Poisson. For this initial data, we obtain a quite differ- 
ent evolution at intermediate times, with different magni- 
tudes of lower-order spherical harmonics generated dur- 
ing the evolution, even though the initial data differs from 
the Kerr-Schild case by a small amount. However, at very 
late times it appears that the scalar field decays accord- 
ing to Eq. (|1.2|) . Our results indicate that the coordinate 
effects discussed by Poisson play an important role in the 
details of the evolution at intermediate times, but they 
do not affect the asymptotic decay rate, presumably be- 
cause of mode mixing in the strong-field region, an effect 
that was not included in Poisson's analysis. 

Because Kerr spacetime is axisymmetric, a 2+1 dimen- 
sional simulation would suffice for the present problem. 
We work in 3+1 dimensions because we have available 
a 3+1 dimensional code (see, e.g. which is 

designed to solve the full nonlinear Einstein evolution 
equations and is being used to study the binary black 
hole problem. This code can be applied to not only the 
Einstein equations, but to any first order strongly hyper- 
bolic system of equations. For such a numerical code it 
is extremely useful to find test problems that are simpler 
than, for example, the binary black hole problem, but 
difficult enough so that they still provide nontrivial tests 
of our numerical algorithms. 

The simulation of late-time tails is just this type of 
problem. Because it is linear and involves fewer dynam- 
ical fields, this problem is simpler than those involving 
dynamical black holes. Yet our treatment of this prob- 
lem contains many of the features currently thought to 
be desirable in a solution of the binary black hole prob- 
lem: wave propagation, multiple computational domains, 
parallelism, black hole excision with no boundary con- 
dition imposed on the excision surface, and constraint- 
preserving boundary conditions on certain fields at the 
outer boundary. These features will be discussed further 
m Section ITTT1 

Although this problem is simpler than the evolution of 
dynamical black holes, it is still technically challenging 
in 3+1 dimensions because of the requirement for high 



resolution, long integration times, and a distant outer 
boundary. As discussed in Section lTTTI we overcome these 
difficulties by the use of multiple computational domains 
and a pseudospectral evolution algorithm. The latter 
yields exponential convergence of spatial numerical er- 
rors for smooth problems, allowing us to achieve a given 
level of accuracy using a small fraction of the compu- 
tational resources that would be required by a unigrid 
finite-difference code. 



II. BASIC EQUATIONS 

A. The Background Spacetime 

We write the spacetime metric in the usual 3+1 form 

ds 2 = -a 2 dt 2 + gij (dx l + dt)(dx j + (3 3 dt), (2.1) 

where is the three-metric, a is the lapse and (3 l is 
the shift. The Klein-Gordon equation will involve these 
quantities and also the extrinsic curvature -Ky, defined 
by 



1 d 



(2.2) 



where is a Lie derivative. 

The Kerr spacetime is expressed in Kerr-Schild coor- 
dinates (t, x, y, z). For a Kerr black hole with angular 
momentum aM in the z direction, the 3+1 decomposi- 
tion of the spacetime in Kerr-Schild coordinates is 



, — <\j + 2111,1 ,. 



B 1 



[i + 2Hi t i t y 

2HIH 1 
~ l + 2Hin v 



(2.3) 
(2.4) 

(2.5) 



Kij = -(l + 2Hl t l t ) / [kljdtH + 2Hl {i d t l j) ] 

-2(l + 2#Z¥)~ 1/2 [By {l^Hl 1 ) 

+2H 2 l t l%d\ k \l j) + HlHil^dkH] , (2.6) 

where H and are given in terms of the black hole's 
mass M and its angular momentum aM by 



H 



I, 



Mr 3 Bh 

r BL + a2 Z 2 



r B LX + ay r B hy - ax 



a 2 ' r| T + a 2 ' r B L 



(2.7) 
(2.8) 



BL ~ u ' BL 

and the Boyer-Lindquist coordinate tbl is defined by 

^2 



x 2 + y 2 



1 BL 



z 



= 1. 



(2.9) 



BL 



Here and in the following, the quantity r without a BL 
subscript refers to the Kerr-Schild radial coordinate de- 
fined by 



2 _ 2 i 2 i 2 

r = x + y + z . 



(2.10) 
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In Kerr-Schild coordinates, the event horizon is located 



at 



(m + Vm 1 



(2-11) 

Notice that the horizon is not spherical in these coordi- 
nates. We typically set a — 0.5; in this case the largest 
coordinate sphere contained within the event horizon has 
a radius of 1.87M, and the smallest sphere that is outside 
the Cauchy horizon has a radius of 0.52M. 

As a consistency check, we compare results for a Kerr 
background with a = to results using a Schwarzschild 
background expressed in Painleve-Gullstrand 0, 0] 
coordinates. In these coordinates, the spatial three- 
metric is flat, leading to a simple representation of the 
Schwarzschild solution: 



9ij 
K i4 



V2A//r 3 % - 3 A /M/2r 3 f i f j 
1, 



[3 k = ^2M/rf k 



(2.12) 

(2.13) 
(2.14) 

(2.15) 



where Sij is the Euclidean metric, r is the areal radial 
coordinate (which for a = is the same as the Kerr-Schild 
and the Boyer-Lindquist radial coordinate), f% — Xi/r is 
the Euclidean unit vector in the radial direction, and M 
is the mass of the black hole. The event horizon is located 
at r = 2M. 



B. Klein-Gordon equation 



teristic fields that propagate with well-determined char- 
acteristic speeds with respect to any two-dimensional sur- 
face, such as a boundary. If the normal to the surface is 
£i, then the characteristic fields are 

u* = n±£%, (2.21) 
«? - '!>- i,i"t",. (2.22) 
ifi = ip. (2.23) 

The fields propagate along null rays (coordinate ve- 
locity v l — —f3 l ± at; 1 ), and the other fields propagate 
along the normal to the spatial hypersurface (coordinate 
velocity v l — —pi 1 ). Note that all characteristic fields 
propagate causally. The decomposition into characteris- 
tic fields is invaluable for the purpose of setting mathe- 
matically consistent boundary conditions. At a boundary 
with normal £j, boundary conditions must be imposed 
only on incoming characteristic fields, that is, those hav- 
ing v l ^i < 0. Boundary conditions must not be imposed 
on other characteristic fields. 

Note that the definition of $j, Eq. (|2.17() . becomes a 
set of constraints, 

dip 
dx l 

that must be satisfied at all times. If d — initially and 
the solution is advanced in time by solving Eqs. I|2.18l - 
I2.20fl exactly, then C, will remain zero for all times, 
as long as the boundary conditions are consistent with 
Cj = 0. However, both numerical truncation errors and 
boundary errors can cause Ci to drift away from zero. 
Therefore, tracking the evolution of d provides a test of 
the accuracy of our simulations. 



(2.24) 



We write the Klein-Gordon equation (|l.lfl in first-order 
form by introducing four new variables: 



n = — ( d ^ & d ^ 

a \ dt dx l 



dip 
dx l 



(2.16) 
(2.17) 



In the background given by equation (|2.1|) . the Klein- 
Gordon equation ll.lj l and the commutivity of partial 
derivatives yield the following system of evolution equa- 
tions 



dip 
dt 
dll 
~~dt 



an, 



-gVQjOtj + aKU, 



n//'<i>,. ; • n ;/ '-T';,<h, 



— = '"'I',.., • IP, "II , lln. 



(2.18) 

(2.19) 
(2.20) 



where indicates differentiation with respect to x b . 

The system (12.1814^7^01 is symmetric hyperbolic, so the 
quantities ip, II and <&i may be decomposed into charac- 



III. NUMERICAL METHOD 

A. Computational Domain 

We solve Eqs. (|2. 18112. 20|l in a 3D spherical shell ex- 
tending from a radius r = r m j n lying between the event 
horizon and the Cauchy horizon to some large radius 
T = ?' max . Because all characteristic fields propagate 
causally, placing the inner boundary inside the event 
horizon means that all characteristic fields are outgoing 
(into the hole) there: > 0. Therefore we impose no 
boundary condition at the inner boundary. Typically we 



choose 



= 1.75M. 



The outer boundary must be placed at a large ra- 
dius because the power-law tails of interest are due to 
backscattering of radiation off the background geometry 
at large distances. If we wish to measure the tail contri- 
bution to the scalar field at time t and radius r, then the 
outer boundary must be placed roughly at r max > r+t/2, 
so that the backscattering responsible for the tail con- 
tribution occurs within our computational domain. Be- 
cause determining the decay rate of tails requires evolu- 
tion to approximately t = 600 M, we typically place our 
outer boundary at r max = 300M . 
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To facilitate multiprocessing, the domain is divided 
into concentric subdomains, each a spherical shell with a 
width of 10M. Each subdomain is evolved independently 
except for boundary conditions, so each subdomain can 
be evolved on a different processor, with interprocessor 
communication occurring only at the boundaries. To im- 
pose boundary conditions at an interdomain boundary, 
we set the time derivative of each incoming characteris- 
tic field equal to the time derivative of the corresponding 
outgoing field from the neighboring subdomain. 

B. Solution Technique 

Our numerical methods are essentially the same as 
those we have applied to the evolution problem in full 
general relativity UA E3 • We use a pseudospec- 

tral technique on each subdomain to evolve Eqs. 1)2.181 - 
I2.2U|) in time. Given a system of partial differential equa- 
tions 

J^/(x, t) = Hf (x, t),df(x, t)/dx*), (3.1) 

where / is a vector of variables, the solution /(x, t) is 
expressed as a time-dependent linear combination of N 
basis functions </>(x): 

N-l 

Mx,t)=£)/ fc (t)0fc(x). (3.2) 

fc=0 

Spatial derivatives are evaluated analytically using the 
known derivatives of the basis functions: 

r) N ~ 1 8 

— f N ( X ,t)=J2Mt)g^M^ (3-3) 

fc=0 

The coefficients fk (t) are chosen so that equation (|3.1|) 
is satisfied exactly at N c collocation points selected from 
the spatial domain. The values of the coefficients are 
obtained by the inverse transform 

N e -1 

fk(t)= M*i,t)<f>k(xi)wi, (3.4) 

i=0 

where Wi are weights specific to the choice of basis func- 
tions and collocation points. One can now transform at 
will, using equations (|3.2|l and (|3.4II , between the spectral 
coefficients /fe(t) and the function values at the colloca- 
tion points /jv(xj,t). The differential equations (|3.1() are 
now rewritten, using equations (|3.2II3.4|) , as a set of ordi- 
nary differential equations for the function values at the 
collocation points, 

^/ JV (x i ,t)=ft({/iv(x J -,*)}) ) (3.5) 
where Qi depends on /jy(xj-,i) for all j. 



Equations l|3.5[) are integrated in time using a fourth 
order Runge-Kutta method. Boundary conditions are in- 
corporated into the right-hand side of Eqs. I|3.5|l using the 
technique of Bj0rhus [l||: if P + is the projection opera- 
tor that annihilates all incoming characteristic fields at a 
boundary, and P~ is 1 — P + , then at each boundary point 
i the differential equation l|3.5|l is modified as follows: 

^f N ( Xi ,t) = p+g i ({f N (x j ,t)}) + p-B i ({f N (xi j ,t)}), 

(3.6) 

where P Bi({/.2v( x 3» *)}) encodes the boundary condi- 
tion placed on the time derivatives of the incoming char- 
acteristic fields. 

For computational subdomains with spherical bound- 
aries, it is natural to use spherical coordinates. We 
choose our basis functions to be Chebyshev polynomi- 
als in radius and spherical harmonics in angles. Although 
our basis functions are based on spherical coordinates, we 
choose our dynamical scalar field variables and our grav- 
itational variables to be the Cartesian components, and 
not the spherical components, of the relevant quantities. 
This allows us to use the same angular basis functions 
for all variables without regard to regularity. 

To eliminate high-frequency numerical instabilities 
that sometimes develop during our simulations, we apply 
a filter to the right-hand sides of Eqs. I|3.5|l before incor- 
porating boundary conditions via the Bj0rhus algorithm. 
The filter consists of simply setting high-frequency spher- 
ical harmonic coefficients to zero. The components that 
are set to zero depend on which equation is being solved: 
if ^max is the index of the highest frequency basis function 
Yi m , then typically the largest I retained in the right- 
hand side of the II equation 1)2.19)1 is 3£ roax /2 — 1, and 
the largest I retained in the right-hand sides of the <E>i 
equations 1)2.20(1 is 3^ ma x/2. This is similar to the "3/2 
rule" commonly used to eliminate nonlinear aliasing er- 
rors 0|- No filtering is done for the ip evolution equa- 
tion 12.18fl . and filtering is not performed on the radial 
basis functions. The degree of filtering necessary to ob- 
tain stability depends on both the background geometry 
and the configuration of subdomains, and is not com- 
pletely understood. For example, in some cases no filter- 
ing is needed, and in others it suffices to zero only modes 
with I = ^ max in the $^ equations (|2.2U|) and modes with 
(■ > dax - 1 in the II equation (|2.19|) . 

C. Outer boundary conditions 

The simplest outer boundary condition is obtained by 
setting the time derivatives of the incoming characteris- 
tic fields u~ , u°, and to zero. While this works well 
for a Schwarzschild background in Painlcvc-Gullstrand 
coordinates, for a Kerr background in Kerr-Schild coor- 
dinates, even for a = 0, this boundary condition produces 
strong violations of the constraint Ci at the outer bound- 
ary, even at t = 0. These constraint violations propagate 
inward and grow, eventually dominating the numerical 
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dt 



= -n/r. 



(3.8) 



This is equivalent to assuming the Sommerfeld condition 
ip = f(t—r)/r, for some unknown function /, at the outer 
boundary. In practice, imposing this boundary condition 
proved less critical than imposing Eq. I|3.7[) . This is be- 
cause, as explained in Sec. IIII Al for studying tails our 
integration time is less than the time it takes light to 
travel from the black hole to the outer boundary and 
back again. Therefore, although it is important that the 
outer boundary condition is well-behaved when there is 
no wave there, the magnitude and nature of the reflec- 
tions produced when a wave passes through the boundary 
is largely irrelevant, because the evolution ends before 
these reflections reach the observation point. 



D. Convergence 



FIG. 2: Constraint violation during a scalar field evolu- 
tion for four radial/angular resolutions (the notation R/lmax 
means that we use R radial basis functions and retain angu- 
lar basis functions up to order ^ max ). Plotted is |[C||, where 
(||C||) 2 = (1/4tt) / C l C l dn. The integration is taken over 
the surface r = 11.75M. For comparison, the value of the 
scalar field at t — 300M at r = 11.75Af is on the order of 
10 -8 . This evolution corresponds to the Y20 case shown in 
Figures |S| and |7| 



solution. Because these constraint violations appear as 
oscillations in the variable ip but do not affect the fields II 
and <i>i , we were able to greatly reduce them by changing 
the boundary condition on ift to 

— = ^n + ^,. (3.7) 

This is the same as Eq. I|2.18[) except that ip t i has been 
replaced by <3?i using the constraint (|2.24|) . Exactly this 
type of boundary condition has been used before 20] in 
the field of numerical relativity, where finding methods 
of constructing boundary conditions that preserve the 
constraints is a topic of active investigation 0, 153, 13, 

mm 

Figure |2] shows the norm of the constraint for several 
different resolutions during an evolution of a scalar field 
in Kerr spacetime, using the boundary condition l|3.7|l . 
The constraint violations decrease rapidly with increas- 
ing resolution. For the higher resolutions, the constraint 
violation is small compared to the magnitude of the scalar 
field. 

Even with the use of Eq. I|3.7I) , reflections (with a small 
constraint- violating contribution) occur when an outgo- 
ing pulse of scalar field reaches the outer boundary. The 
reflected pulse then propagates (causally) inward. These 
reflections can be reduced by modifying the boundary 



The convergence properties of a pseudospectral code 
is more difficult to analyze than, for instance, a second- 
order finite-difference code. This is because in the former 
there are several sources of truncation error that scale dif- 
ferently with resolution. Spatial truncation errors should 
converge exponentially (and errors associated with the 
radial direction may have a different exponential conver- 
gence rate than those associated with the angular direc- 
tions because different basis functions are used). Time 
integration errors should scale like (At) 4 because we are 
using a fourth-order Runge-Kutta time integrator. Fur- 
thermore, the Courant condition constrains At as a func- 
tion of spatial resolution. For the resolutions we use, the 
scaling is roughly max(Ai) ~ N~ 2 , where iV r is the num- 
ber of radial collocation points; the scaling is not simply 
At ~ NjT 1 because the collocation points are distributed 
nonuniformly. 

Figure[3]shows the norm of the constraint as a function 
of radial resolution at different times, for the evolutions 
shown in Figure [3 The angular resolution is fixed but 
At is varied so that the Courant condition remains sat- 
isfied. The convergence is exponential, indicating that 
the radial spatial errors dominate both the time integra- 
tion errors and the angular integration errors. Even at 
late times, when the scalar field is very small, the con- 
vergence plot is still roughly exponential, although it is 
noisier than at early times. Figure 0] shows the norm of 
the constraint as a function of angular resolution at differ- 
ent times, for fixed (high) radial resolution and fixed At. 
At late times, the convergence is exponential for low res- 
olution and then saturates when the angular truncation 
error drops beneath radial truncation error. For early 
times the angular truncation error is already small, even 
at low resolution (as is expected for initial data that is 
pure 1 = 1 without higher angular components), so ex- 
cept for the difference in the two lowest resolutions one 
does not see any dependence on the number of angular 
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FIG. 3: The same constraint norm as shown in Figure 
but plotted versus the number of radial collocation points for 
three different values of t. The angular resolution is fixed at 
4u = 17. These evolutions correspond to the 120 case shown 
in Figures [S] and |7| 



FIG. 4: The same constraint norm as shown in Figure [5] but 
plotted versus the angular resolution £ max for three different 
values of t. The radial resolution is fixed at N r = 36 and At 
is fixed at 0.055M. These evolutions correspond to the Y20 
case shown in Figures El and Q 



collocation points. Presumably, for large enough spatial 
resolution, the fourth-order time integration error should 
dominate (unless the errors drop beneath roundoff level 
first), but we do not see this for the resolutions considered 
here. 



IV. RESULTS AND DISCUSSION 

A. Schwarzschild 

As a test of our numerical techniques, we began by 
evolving the well-understood case of a scalar field in a 
Schwarzschild background. We write the background in 
Painleve-Gullstrand coordinates, and we choose initial 
data of the form 



$ = 0, 
*i - 0, 
n = n o (r,0,< 



e -(r-ro)V" 2 y, oroo (0,<£). 



(4.1) 
(4.2) 

(4.3) 



Figure n shows the results of a simulation with £q = I, 
too = 0. For all simulations shown here we set ro = 
12M and w = 2M. Plotted are results obtained using 
resolution 24/9 (where the notation R/£ max means that 
we use R radial basis functions and retain angular basis 
functions up to order £ max )- From roughly t = 40M to 
140M the scalar field behaves like ip ~ e~ iuJt with ujM ~ 
0.29 — 0.097i. This agrees with published values 26] of 
the least-damped quasinormal frequency for scalar I = 1 
perturbations of the Schwarzschild geometry. 



As time increases, the decay of the scalar tails ap- 
proaches the expected power-law decay ip oc . Since 
we cannot numerically evolve the scalar field out to in- 
finite time, our results do not exactly match the analyt- 
ically predicted power law. To facilitate the determina- 
tion of the power law governing the decay of the scalar 
field during the tail phase, the scalar field and its time 
derivative were combined into a "numerical power index" , 
following : 



-t 



UN 



L'l 



where the L2 norms are defined by 



L2 ) 2 = (i/4^) / fdn, 



(4.4) 



(4.5) 



with integration over a surface of fixed r. We compute 
ip using Eq. H2.16|l rather than taking numerical time 
derivatives of ip. 

At finite times during the tail phase, the scalar field 
behaves like 



iP(xt-^+0(t-^ +1 ^) + ... 
which implies that 



(4.6) 



(4.7) 



At late times, the power index asymptotically approaches 
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t/M " t/M 



FIG. 5: The top graph shows the evolution of the power index 
Hn evaluated at r = 11. 9M as a function of time for the case 
shown in Figure Q and for the same case run at a lower res- 
olution. The bottom graph shows differences in /xjv between 
runs done at different resolutions. 

Figure [5] shows the evolution of the power index for 
two different radial resolutions. The power index is ap- 
proaching a value of five, which corresponds to the pre- 
dicted t~( 2 '+ 3 ) decay rate. Moreover, as the resolution 
increases, the numerical results converge, as can be seen 
from the bottom graph in Figure |SJ 

According to Eq. I|4.7|l , we can obtain a better estimate 
for /i by performing a least-squares polynomial fit to \x n 
as a function of t~ l . The least-squares fit will also give 
us an error estimate 28] for /i as long as we provide error 
estimates for each of our numerical values of hn- Our 
code provides fiN at a discrete set of time values fj. To 
estimate the error in /i/v at time ti we use 

Shn(U) = max {\n N {tj\ hi) - (j, N (tj; lo)|} , (4.8) 
je{\i-j\<W\ 

where HN{tj] hi) and /iN(tj]lo) represent numerical val- 
ues for ii n at the highest and next-highest spectral res- 
olution that we used. The purpose of maximizing the 
error over neighboring points is to treat the cases in 
which /ipjitj; hi) and i^N(tj', lo) spuriously agree at a sin- 
gle point — without the maximization this point would 
have an artificially small error estimate. The size W of 
the maximization window is typically 1/20, where / is 
the total number of discrete values of /in(U) that we use 
for the fit. Near the points to and ti we translate the 
maximization window in Eq. Q4.8p . so that, for example, 
for i = the window goes from j = to j = 2W + 1. 

Using a linear fit to the form (|4.7|l we obtain /i = 4.99± 
0.01, and using a quadratic fit we obtain fi = 5.00 ±0.08. 
These agree with the accepted value to within about a 



FIG. 6: Evolution of a scalar field for three different cases: 
initial data proportional to Ybo, Y20, and Y40. Plotted are the 
L2 norms of ip evaluated on the surface r = 11.75M. The 
resolution is 20/12 for £ = 0, 24/17 for £ = 2, and 27/26 for 
£0 = 4. Higher resolution is needed for larger £0 in order to 
resolve the much smaller late-time tail. 

percent. We perform the fits only for data in the tail 
region of Figure GO that is, for t > 400M. The estimate 
of /1 is relatively insensitive to the exact region of t in 
Figure [5] that we choose to perform the fit. 

B. Kerr 

Following our numerical trials with a Schwarzschild 
background, we turned our attention to scalar fields 
around rotating black holes. For our background space- 
time we used a Kerr geometry with spin a — Q.bM . Fig- 
urc|Hldisplays the evolution of the scalar field on this Kerr 
background for three different choices of initial data. The 
initial data are taken to have the form l|4.3|l . where now 
(r, 8, 4>) are related to the Kerr-Schild coordinates (x, y, z) 
defined in Section III Al in the usual way 

x = r sin cos 0, (4-9) 
y = r sin 6 sin (f>, (4-10) 
z = rcosd. (4-11) 

The top, middle, and bottom plots in the figure show the 
decay of a scalar field initially proportional to loo 7 Y20, 
and I40, respectively. Since the latter two cases have 
the same initial value of mo = and are of even par- 
ity, we expect that a Yqq mode will be generated during 
these evolutions. Thus, according to the simple mode- 
mixing picture discussed above, all three evolutions in 
Figure El should follow the same power-law decay at late 
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FIG. 7: The evolution of the power index fiN evaluated at 
r = 11.75Af for the same cases as Fig. E] The power indices 
appear to asymptotically approach a value of three. For each 
case, two resolutions are shown to demonstrate convergence. 
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5.23 ±0.19 
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KS 


0.5 


3.001±0.003 
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BL 


0.5 


2.8 ±0.3 



TABLE I: Numerically-determined power-law decay rates. 
Shown are the spherical harmonic indices of the initial data, 
the (r, 6*) coordinates used in the initial data function Pio 
in Eq. 14.31 (Painleve-Gullstrand, Kerr-Schild, or Boyer- 
Lindquist), the spin of the hole, and the best-fit power index 
at late times. 

For example, initial data proportional to I21 is forbidden 
to evolve by mode-mixing into any lower Yt m mode {£ = 
is forbidden by to conservation and £ = 1 is forbidden by 
parity). Therefore, the i~( 2i + 3 ) law predicts a tail decay 
rate of t , which is what we observe. 

C. Coordinate effects 



times. Figure El supports this prediction: although the 
quasinormal ringing phases are dissimilar, the slope of 
the tails do appear to match. 

The evolution of the power index for these cases is 
shown in Figured The power indices approach a value of 
three, which corresponds to the Price decay rate formula 
ip ~ t~( 2£ + 3 ) for an £ — mode. Estimates of fj, obtained 
by least-squares fits to the numerical data can be found 
in Tabled These estimates all fall within a few percent 
of the value fi = 3. 

The tails have a much smaller magnitude for evolutions 
with larger initial values £$. For example, for £q = 4, 
the scalar field approaches a magnitude of 10~ 12 at late 
times, and its time derivative is two orders of magnitude 
smaller. Thus we are forced to use larger number of spec- 
tral coefficients to resolve the large £q cases. However, be- 
cause the accuracy of pseudospectral methods increases 
exponentially with the number of collocation points, in- 
creasing the number of coefficients only by roughly a fac- 
tor of two in each dimension was sufficient to resolve even 
the £o = 4 case. 

Note that the £0 — 4 case shown in Figures and [7| 
is the smallest value of £q for which the two analytical 
predictions, Eqs. I|1.2fl and i|1.3|) . disagree. Our results 
support the simple picture leading to (|1.2[) . which yields 
a t~ 3 falloff for this case, rather than Eq. I|1.3|) . which 
predicts a t~ 5 falloff. 

Table [I] summarizes various cases that we have studied 
numerically. In addition to the too = cases discussed 
above, we have also computed power-law decay rates for 
too = 1 and too = 2. These cases allow us to test the 
predictions of Eq. 1|1.2|) and Eq. (|1.3fl more thoroughly. 



It has been argued 0] for a rotating weakly-curved 
spacetime that the difference between the predictions of 
Eq. (|1.2|) and Eq. I|1.3f) is related to the choice of coor- 
dinates, and in particular, that Eq. I|1.3|) is correct if the 
initial data were proportional to a spherical harmonic in 
spheroidal coordinates. To test whether this might be 
true for Kerr, we repeated the evolution of initial data 
proportional to I40, but with a small coordinate change: 
We still take initial data of the form (|4.3|) , but we choose 

n = n o (r B L,0BL,</>), (4.12) 

where #bl is the Boyer-Lindquist coordinate defined by 
cos^bl = -z/reL- Because £?bl 7^ 9 and tbl 7^ r, this 
initial data differs slightly from the form (|4.3|) . In fact, 
the magnitude of this difference is only about one part 
in a thousand. 

For brevity, here and in the following we will refer to 
the evolution of initial data (|4.12|) as the "BL" case, and 
we will refer to the £0 = 4 evolution shown in Figure 
as the "KS" case, since the two evolutions differ only 
in which radial and polar coordinate (Boyer-Lindquist 
or Kerr-Schild) are used in the expression (|4.3|l for the 
initial data 3 . Note that the Kerr background is expressed 



3 In principle, instead of evolving this BL initial data, we could 
have expanded the data in terms of KS spherical harmonics and 
suitable radial basis functions, and used linearity to compute 
the result. However, this would require knowledge of both the 
power-law decay rate of each KS spherical harmonic and the 
mixing rates between all pairs of KS spherical harmonics. 
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FIG. 8: Scalar field evolution in a Kerr (a/M = 0.5) back- 
ground resulting from initial data proportional to Y4o(6bl, <j>) 
(labeled by "BL"), plus a corresponding evolution resulting 
from initial data proportional to Kto(#, <f>) (labeled by "KS"). 
Shown are the absolute values of selected Kerr-Schild spectral 
coefficients of ip for both cases. For both cases the resolution 
is 30/29. The initial data for the KS and BL differ by only 
0.1%, yet the details of the evolutions are quite different. 



in the same coordinate system for both the BL and the 
KS evolutions, and is given by Eqs. l|2.3H2.t)[) . Note also 
that the transformation between the BL and KS radial 
and polar coordinates is exactly the transformation (|1.4h 
[T3jl studied in @. 

If the argument of 0] applies to the Kerr geometry, 
then for the BL case, the power law falloff rate should be 
i~ 5 , in agreement with Eq. I|1.3|) . rather than t~ 3 , which 
is predicted by Eq. (|1.2[) . 

It is quite difficult in the BL case to obtain an accu- 
rate value for the late-time power index /u, by the method 
used in sections IIV Al and IIVBI This is because even 
though the initial data l|4.12[l differs only slightly from 
that used in the KS case, the evolution proceeds quite 
differently: The resulting late-time tail is a few orders of 
magnitude smaller than for the KS case, and by the time 
the solution displays its late-time asymptotic behavior, 
the scalar field time derivative is so small (~ 10 -16 ) that 
machine roundoff error (not numerical truncation error) 
obscures the results. 

Fortunately, this roundoff error turns out to be largest 
at high angular frequencies, so it is still possible to de- 
termine the late-time behavior of the BL case for low- 
frequency spherical harmonic components of the solution. 
In Figures |5] and ^\ different spherical harmonic compo- 
nents of the solution are plotted as a function of time for 
both the BL and KS evolutions. The spherical harmonic 



FIG. 9: Same as Fig. |H| showing detail at late times. 

components are computed by 

il> tm = J i>Y em (9,(t>)<m, (4.13) 

where the integral is taken over the surface r = 11.75A/ . 
Note that for all plots shown in Figures |H1 and the 
spherical harmonic appearing in the integral H4.13(l is de- 
fined using the Kerr-Schild 9 and <j> coordinates, and the 
integral is taken over a surface of constant r, not a surface 
of constant tbl- Thus the quantities plotted in Figures |HI 
and are in all cases the spherical harmonic components 
of if> with respect to the Kerr-Schild coordinates. Note 
also that ip£ m and the analogous quantity ipt m can be 
used to compute a power index for each individual spher- 
ical harmonic component. 

For the KS evolution shown in Figures |H1 and [51 the 
initial data consist of pure I = 4, but I = 2 and £ = 
components appear at very early times because of mode 
mixing. The tail of the £ = component can be seen as 
early as t = 150M, but does not exceed the quasinormal 
ringing of the £ — 4 component until t = 300M, after 
which it dominates. Its measured power index is fi — 
3.001±0.003. The tail of the I = 2 component can be seen 
for t > 300Af, but it is extremely small (~ 10~ 14 ). Its 
decay rate is roughly t~ 7 , but it is difficult to determine 
the exponent accurately because it is buried in the noise. 
The tail of the £ = 4 component cannot be seen; the 
£ = 4 component is buried in machine roundoff error 
after t = 400M. 

The BL evolution shown in Figures 00 and |^1 is initially 
almost identical to the KS case. Initially the BL case is 
not pure I = 4 (recall £ here refers to the index of the 
Kerr-Schild harmonic; the BL case is pure t = 4 with 
respect to Boyer-Lindquist spherical harmonics), but also 
has a very small mixture of other components, the largest 
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being I = 6 (not shown) and 1 — 2. As the BL case 
evolves in time, the tail of the I = component first 
appears at t = 200 M, but it is a few orders of magnitude 
smaller than the corresponding I = tail for the KS case. 
Its power index is fi = 2.8 ± 0.3. The tail of the £ = 2 
component, however, appears at t — 250M and is a few 
orders of magnitude larger than the corresponding £ = 2 
tail for the KS case. Its power index is [i = 7.0 ± 0.5. As 
in the KS case, the tail of the £ = 4 component cannot 
be seen because of machine roundoff error. 

Although at intermediate times the £ = 2 mode is im- 
portant for the BL case, it is clear from Figure [5] that 
at very late times, the £ = mode will eventually dom- 
inate, resulting in a decay rate of t~ 3 . In other words, 
the asymptotic decay rate appears to be independent of 
whether the initial data are expressed in terms of Kerr- 
Schild or Boyer-Lindquist spherical harmonics. Thus, the 
argument in Q apparently does not carry over to the 
Kerr geometry. This is presumably because Kerr has 
strong-gravity regions that influence the scalar field, and 
strong-gravity effects were not included in • 

We perform all our evolutions of BL initial data us- 
ing a background expressed in KS coordinates; this is 
because black hole excision requires coordinates that are 
regular through the horizon. A natural question to ask 
is whether our results are different than they would be 
if we had expressed our background in BL coordinates. 
The answer is yes, but only because we set "BL initial 
data" on a hypersurface of constant KS time, not on a 
hypersurface of constant BL time. Setting initial data 
on a hypersurface of constant BL time would require an 
integration in BL coordinates (at least until the solution 
were known on a full hypersurface of constant KS time, 
and from that point on the evolution could proceed in 
KS coordinates), and is beyond the scope of this paper. 
However, for investigating the coordinate effect described 



m , it is unnecessary to evolve initial data on a BL time 
slice; the derivation of this effect in (2| involves no change 
in time slicing but only a transformation of spatial coor- 
dinates, the same transformation that we have done here. 

It would also be interesting to repeat the evolutions 
in this section with an outgoing initial pulse centered far 
from the black hole, so that only the weak-gravity re- 
gion is seen by the scalar field until extremely late times. 
In this case the weak-field approximation assumed by 
would be valid for an extended period of time, so dur- 
ing this time one should sec a difference in decay rates 
between the BL and KS evolutions. Such a computation 
would be more difficult than the ones presented here be- 
cause it would require a more distant outer boundary and 
therefore longer integration times. A single run such as 
the one shown in Figures 0D and [§] takes about 23 hours 
on 30 processors of the IA-32 Linux cluster at NCSA; 
the run time scales like if all subdomains have the 
same number of radial collocation points N r . Future 
work may involve a self-gravitating scalar field; in this 
case the equations would be fully nonlinear. 
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